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ABSTRACT 


A two-dimensional numerical model is developed for the 
unsteady oscillatory combustion of the solid propellant 
flame zone. Variations of pressure with low and high 
frequency responses across the long flame, such as in the 
double-base propellants, are accommodated. The formulation 
is based on a premixed, laminar flame with a one-step 
overall chemical reaction and the Arrhenius law of 
decomposition for the gaseous phase with no condensed phase 
reaction. Numerical calculations are carried out using the 
Galerkin finite elements, with perturbations expanded to 
the zeroth, first, and second orders. The numerical 
results indicate that amplification of oscillatory motions 
does indeed prevail in high frequency regions. For the 
second order system, the trend is similar to the first 
order system for low frequencies, but instabilities may 
appear at frequencies lower than those of the first order 
system. The most significant effect of the second order 
system is that the admittance is extremely oscillatory 
between moderately high frequency ranges. 


1 



TABLE OF CONTENTS 


Page 

ABSTRACT i 

TABLE OF CONTENTS ii 

LIST OF FIGURES iv 

LIST OF TABLES ..;.... vii 

LIST OF SYMBOLS. .-. . . viii 

I. INTRODUCTION 1 

II. PRELIMINARIES FOR COMBUSTION WAVES 3 

2.1 Prexnixed Laminar Flame 3 

2.2 Fluid Dynamically Excited Oscillations.... 6 

2.3 Propellant Characteristics 6 

2.4 Acoustic Admittance 7 

III. COMBUSTION DYNAMICS 10 

3.1 General 10 

3.2 Governing Equations 10 

3.2.1 Gas Phase 11 

3.2.2 Solid Phase 13 

3.2.3 Perturbation Expansions 15 

3.3 Finite Element Applications 27 

IV. DISCUSSIONS 31 

V. CONCLUSIONS 43 

REFERENCES 45 

ii 


APPENDICES 


A. DERIVATION OF NON-DIMENSIONAL GOVERNING 

EQUATIONS 79 

B. PERTURBATION OF MODEL EQUATIONS 88 

C. DERIVATIONS OF BOUNDARY EQUATIONS 93 

D. DERIVATIONS OF FINITE ELEMENT EQUATIONS 112 


ill 



LIST OF FIGURES 


Figure Page 

3.1 Domain of study for one-dimensional case... 47 

4.1 Finite element geometry 48 

4.2 Steady-state distributions of field 

variables in flame zone 49 

4 . 3 First order pressure distributions 

vs. frequency 50 

4.4a First order distributions of field 

variables at « = 1 (real parts) 51 

4.4b First order distributions of field 

variables at « * 1 (imaginary parts) 52 

4.5 First order density distributions 

vs. frequency 53 

4.6 First order temperature distributions 

vs. frequency 54 

4.7 First order species (fuel) distributions 

vs. frequency 55 

4.8 First order velocity distributions 

vs. frequency 56 

4.9a First order acoustic admittance and 

burning rate vs. frequency 57 


IV 



4.9b Steady oscillatory mode: typical dependence 

of amplitude on frequency of pressure 
oscillation 58 

4.9c Real parts of acoustic admittance and 

burning rate vs. frequency 58 

4.10 Second order time-independent pressure 

distributions vs. frequency 59 

4.11 Second order time-independent distributions 

of field variables at u = 1 60 

4.12 Second order time-independent density 

distributions vs. frequency 61 

4.13 Second order time-independent temperature 

distributions vs. frequency 62 

4.14 Second order time-independent species 

(fuel) distributions vs. frequency 63 

4.15 Second order time- independent velocity 

distributions vs. frequency 64 

4.16 Second order time-dependent pressure 

distributions vs. frequency 65 

4 . 17a Second order time-dependent distributions 

of field variables at « ■ 1 (real parts) ... 66 

4 . 17b Second order time-dependent distributions 
of field variables at u = 1 (imaginary 
parts) 67 

4 . 18 Second order time-dependent density 

distributions vs. frequency 68 


v 



4 . 19 Second order time-dependent temperature 

distributions vs. frequency 69 

4.20 Second order time-dependent species 

(fuel) distributions vs. frequency 70 

4.21 Second order time-dependent velocity 

distributions vs. frequency 71 

4.22 Second order pressure distributions 

vs. frequency 72 

4.23 Second-order distributions of field 

variables at u - 1 73 

4.24 Second order density distributions 

vs. frequency 74 

4.25 Second order temperature distributions 

vs. frequency 75 

4.26 Second order species (fuel) distributions 

vs. frequency 76 

4.27 Second order velocity distributions 

vs. frequency 77 

4.28 Offset burning rate and acoustic admittances 

of the second order time-dependent and 
time-independent systems 78 


vi 


LIST OF TABLES 


Table Page 

1 Typical value and range of parameters 

used in numerical calculation 32 


vii 


LIST OF SYMBOLS 


a speed of sound 

B frequency factor for gas phase reaction 
C p specific heat at constant pressure for gas 
C 8 specific heat of solid 

Da Damkohler number 

E gas phase activation energy, acoustic energy in 
Eq. (2.14) 

E s surface activation energy 

F eigenvector 

G source term in numerical formulation 
h heat of combustion per unit mass 

H heat of reaction 

I order of perturbation 

k thermal conductivity, constant in Eq. (3.11), 

wave number (2.29) 
k reaction rate constant 

ft characteristic flame length 

L latent heat of vaporization, chamber length in Eqs. 
(2.25) 

m mass flux 

M b Mach number 

n order of chemical reaction 

P pressure 


vm 



Pr Prandtl number 

Q dimensionless heat of gas phase reaction 
r burning rate 

R gas constant 

R c dimensionless distance from surface in Eq. (3.22) 

t time 

T temperature 

u axial gas velocity 

Uj velocity vector 

u c core velocity in Eq. (3.22) 

v normal gas velocity 

V total volume 

w reaction rate 

W molecular weight 

x axial distance parallel to the surface 
X solution vector, mole fraction 

y normal distance from the surface 

Y mass fraction of fuel species 

z oxidizer-fuel ratio, z-axis in cylindrical coordinate 
a thermal diffusivity 

fl ratio of solid to gas density 
T ratio of heat capacity, surface area 

5 temperature exponent in Eq. (3.7) 

7 specific heat ratio 

e perturbation parameter 

f = Jc 8 *C p */k*C g * 


xx 


X 


eigenvalue 
X 7 Lagrange multiplier 
u viscosity coefficient 
v stoichiometric coefficient 

% - k*/k,* 

p density 

<t> scalar potential in Eq. (2.26), interpolation 

function 
« frequency 

ft interior volume 

Subscripts and Superscripts 
* dimensional quantity 

- steady-state mean value 

A spatial variable 

(i) i th order perturbation coefficient 

+ gas side of interface 

solid side of interface 
c propellant cold side 

D time-dependent variable 

f fuel 

h homogeneous solution 

i,j vector quantity 

I time-independent variable 

0 steady-state, oxidizer 

p particular solution, product 

r reference value 


x 




s 


solid phase 


a. , & , 


7, a finite element global node number 
t complex conjugate 


xi 


CHAPTER I 


INTRODUCTION 

Determination of oscillatory pressure and velocity 
variations from their mean values in solid propellant 
combustion chambers has been the subject of study for many 
years. Because of computational difficulties, however, 
rigorous and accurate methods of solution still remain in 
various stages of development. Often details of physics 
are obscured in complicated computational strategies. In 
view of the fact that reliable experimental measurements 
are difficult to obtain, it is crucial to possess an 
analytical tool to accurately predict combustion dynamics 
and combustion efficiency in designing successful solid 
propellant rocket motors. 

Earlier works on combustion dynamics include Hart and 
McClure [1], Denison and Baum [2], Culick [3] and T'ien [4] 
among others. These studies are centered around one- 
dimensional, linear oscillatory burning. Culick [5,6], 

Baum and Levine [7], and Yang et al. [8] subsequently 
studied the nonlinear growth and limiting amplitude of 
acoustic waves in a combustion chamber. Nonlinear behavior 
as characterized by higher order perturbation expansions of 
governing equations has been investigated by Flandro [9] 
and Chung and Kim [10] in recent years. They assumed that 
the pressure varies sinusoidally with time, the pressure 
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amplitudes being a function of position. 

In the present approach, the pressure is assumed to 
vary in the unstable state. This allows the presence of a 
long flame such as in double-base solid propellants, in 
which high pressures and high frequencies may be 
accommodated. Galerkin finite elements [11] are utilized 
to model numerically the geometries of burning surfaces and 
the flame thickness. The nonsteady governing equations are 
linearized by means of the zeroth, first, and second order 
perturbation expansions. The boundary conditions, 
including the burning surfaces and flame edges, are imposed 
by means of the Lagrange multipliers. 

In the simple example problems demonstrated in this 
paper, the gaseous flame is assumed to follow the Arrhenius 
law with one-step forward chemical reactions. No condensed 
phase reaction is included in the formulation. The 
calculated results confirm the prediction reported by other 
investigators in the literature for the low frequency 
region. Extended studies are then carried out for the high 
frequency region. It has been found that amplification of 
oscillatory motions does prevail in the high frequency 
region. Discussions are presented in Section 4. 

The present study does not include turbulent boundary 
layers which may play an important role in erosive burning 
and combustion instability; this is the subject of the 
current study. 



CHAPTER II 


PRELIMINARIES FOR COMBUSTION WAVES 
2.1 Premixed Laminar Flame 

Combustion is the chemical process which gives rise to 
the conversion of reactants into products. During 
combustion, chemical reactions are coupled with heat and 
mass transfer. Energy and species are transported with 
rapid chemical reactions occurring exothermically in single 
or composite phase. Most chemical reactions occur in a 
chain sequence composed of four basic processes: 
initiation, branching, propagation, and termination of the 
radicals. A chemical reaction to combustion can be written 
in the form 


N 

E 

a=l 


V » 

a , to 


N 

l r" Y a 

a=l a ' m 


( 2 . 1 ) 


where Y a denotes the chemical species a and y a ' /Tn and y a " /In 
represent the stoichiometric coefficients of species Y a for 
the reactants and for the reaction products in the m* h 
reaction, respectively. The stoichiometric coefficients 
v a ' t m and r a ” m must satisfy elemental material balances. 

The phenomenological principle of mass action states that 
the rate of reaction is proportional to the product of the 
concentration of the reactants. Thus, the reaction rate is 
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given by 


w 


a 


-W 


a 


M 

E 

m=l 


(y" 
a , m 


y • 
a / m 


>*■ 


N 

n 

b«i 


'V' 

iRT > 


v • 

6 , m 


(2.2) 


in which is the mole fraction and k m is the reaction 
rate constant, with m indicating the reaction order. The 
order of a reaction is determined by summing the exponents 
of the reactant mass fraction that express the reaction 
rate. The constant k m for a single step reaction is 
empirically given by the Arrhenius law 


“■m 


A m exp 


V 


RT 


( 2 . 3 ) 


with R being the universal gas constant, E ra the activation 
energy, and A m the frequency factor for the reaction step. 
According to the theory related to the collison process and 
potential energy barrier of the reaction, E m is the minimum 
barrier height which reactants must acquire before reaction 
process. The factor exp( _ E 1 „/RT) represents the fraction of 
collisions between reactant molecules in which products can 
be formed. Generally, E m is constant, but A m depends on T 
such that 


A. - B.T ! * 

where B m and a m are constant. Combining those relations 
gives 
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w 


a 


-W 


a 


M 

I 

m=l 


(/•' 
a / m 


v * )B,.T*ine' 
a , rn 


N 

E m/ RT > n 

11=1 


'Vj 

iRT J 

(2.5) 


in which the detailed reaction kinetics are neglected. 

For a single step irreversible chemical reaction, the 
chemical process may be written in the form 


y'F + k'O 4 y"P (2.6) 

f o p . 

with F, 0, and P -being the species representing the fuel, 
oxidizer, and product. In this case, the overall reaction 
rate is given by 


dY f ( E 

Wf = = -BT s Y f r Y 0 s exp 

dt RTj 


(2.7) 


It is noted that the exponents r and s are not necessarily 
related to the stoichiometric coefficients y f ' and y 0 '. In 
many cases, actual values of B, 5, r, s, and E are 
empirically determined. Based on this model, the local 
heat release rate Q f can be determined from Eq. (2.7) and 
the heat of combustion per unit mass of fuel h f , 


Qf 


h f dY f 
P dt 


( 2 . 8 ) 


where Q f represents heat addition to the system. Thus, 
Eqs. (2.7) and (2.8) can be used in conjunction with the 
species and energy transport equations, respectively. 
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2.2 Fluid Dynamically Excited Oscillations 

The contribution of flow excitation is not usually 
distinguishable from that due to combustion. However, 
there are three rather distinct examples that indicate flow 
excitation by means other than combustion: (1) oscillations 
related to vortex shedding and transport from the segment- 
transition regions of the motor cavity in large segmented 
motors at low frequencies [13]; (2) excitation by the 
interaction between transverse flow oscillation and axial 
mean flow [14]; and (3) excitation by turbulence [15]. The 
contribution of flow excitation to oscillations is usually 
small compared to the combustion oscillations, primarily 
because flow excitation results from the kinetic energy 
flux in the mean flow, while combustion oscillation results 
from the much larger thermal energy flux in the combustion 
zone. However, flow excitation may play a role in 
initiation of combustion-driven oscillations and 
correspondingly affect stability limits. 

2.3 Propellant Characteristics 

Experience shows that combustion instability is highly 
dependent on propellant characteristics. It is clear that 
the effects can be due to either changes in combustion 
dynamics of the propellant or changes in damping. Although 
the stability is a property of the entire combustor and 
cannot generally be attributed to the propellant alone, the 
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variables involved in the response function are related to 
the characteristics of propellants. The growth constant in 
the response function shows that a high energy, fast 
burning propellant in a charge configuration with a large 
burning surface area will tend to be unstable unless the 
real part of the response function is small or the damping 
effect is large. It has been shown that the value of the 
response function for energetic propellants is high enough 
at some frequencies to cause motor instabilities. These 
propellants include ammonium perchlorate, double-base 
propellant, aluminized propellants, and combustion of 
propellants with other energetic oxidizers. 

2 . 4 Acoustic Admittance 

The acoustic wave equation for an adiabatic 
compressible inviscid flow is given as 

1 3 2 t 

0 (2.9) 

a 2 3t 

where a is the speed of sound and ^ is a scalar potential 
defined as Uj = 

The corresponding acoustic pressure P* is 
3 4 > 

p. = - p — (2.10) 

3t 

Taking a partial derivative of Eq. (2.10) leads to 
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dU { 1 

- P' (2.11) 

3t p ' 1 

The general solution of Eq. (2.9) in the one-dimensional 
case is 

^ = Ae i<kx - Wt> + Be- i(kx+ut) (2.12) 

in which A and Bare complex quantities, k is the wave 
number, and u is-the frequency. If the positive going 
travelling wave is considered, it follows that 

u = ce i,kx " wt) (2.13) 

Substituting Eq. (2.13) into Eq. (2.11) yields 

P' = Cpae* { kx_wt » (2.14) 

From Eqs. (2.13) and (2.14), we have 

P' 

— -pa (2.15) 

u 

The quantity pa is called the specific impedance in 
acoustics and the reciprocal of the normalized specific 
impedance is denoted as the admittance in combustion 
instability. Thus, the acoustic admittance Y' is given as 

u 


Y' 


P* 


(2.16a) 



Dividing by the steady-state value gives the non- 
dimensional form of y' ( 

P 0 u 

Y = (2.16b) 

u 0 P* 

where the subscript o indicates the mean flow. 

The so-called response function F is defined as 

F - ReCm'/moJ/^P'/P'o) (2.17) 

where m is the mass flow rate and Re is the real part of 
the complex quantity. Both acoustic admittances and 
response functions are representative of the combustion 
instability. In what follows, however, we shall use the 
acoustic admittance as a measure of combustion wave 


instabilities. 



CHAPTER III 


COMBUSTION DYNAMICS 


3. 1 General 

A numerical model is developed for the unsteady state 
combustion of a solid propellant subject to acoustic 
pressure oscillations in both low and high frequency 
ranges. The governing equations are perturbed to the 
second order and then solved by using the finite element 
method. After reproducing the analytical model given by 
Friedly and Petersen [16,17], the numerical results are 
verified for the one-dimensional case by comparison with 
the analytical solution. Thus, the work is extended for 
the multi-dimensional case where the critical assumptions 
required in the analytical method may be removed. As a 
result, the investigation of variable distributions in the 
flame zone is made possible and, consequently, the acoustic 
admittance is obtained to determine the combustion 
instability in a more realistic manner. 

Discussions of governing equations and finite element 
solutions will be presented in the following sections. 

3 . 2 Governing Equations 

The governing equations for a premixed laminar flame 
are derived under the following assumptions: (l) the Mach 
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number is small, (2) body forces are negligible, 

(3) radiative heat transfer is negligible, (4) the Soret 
and Dufour effects are negligible, (5) diffusion caused by 
the pressure gradient is negligible, (6) the gas mixture is 
ideal and thermal coefficients are constant, and (7) the 
Lewis number is unity for the explicit relation between 
concentration and temperature. A schematic representation 
of the coordinate system for a one-dimensional model is 
shown in Fig. 3.1-. 

3.2.1 Gas Phase 

For the simplicity of the combustion modeling of solid 
propellants, the gaseous flame is assumed to be multi- 
dimensional, premixed laminar, and calorically perfect, and 
a one-step forward chemical reaction occurs. The 
combustion of a solid propellant is approximated by the 
Arrhenius law. Thus, the conservation laws for the multi- 
component reactive system in the gas-phase are represented 
in the non-dimensional forms as follows (see Appendix A for 
details) : 

Continuity 

dp 

— + ( P u i ) , i ■ 0 


(3.1) 
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Momentum 


9u i 

+ pu i/jUj 



Pr 


u 


i / j j 


1 ' 
+ “ U j,ij 


Energy 


0 

(3.2) 


9T 7 - 1 3P 

p — + pU i T y i T /U - wh = 0 (3.3) 

at y 3t 

Species 


ay 

P — + pU i Y / j 

at 


y /U + W = 0 


(3.4) 


State 

P ■ pT (3.5) 

where the commas denote partial derivatives, the repeated 
indices imply summing, Pr is the Prandtl number, and Y 
represents the fuel mass fraction. Note that only one out 
of three species equations (fuel, oxidizer, and product) is 
taken into account from the simple chemical reaction model 
[4]. The following characteristic parameters are used to 
render the above equations dimensionless: 

P = p*/ Po * , P = P*/P 0 * / T = T*/T 0 * 

Ui - u i */v 0 * , t = t*v 0 */j!* / Xi = Xj */jt* (3.6) 
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M b - v 0 */a 0 * # h = h*/C p *T 0 * , w = w*a*/ p 0 *v 0 * 2 

in which jt* is the flame thickness given by a*/v 0 * , with a* 
being the thermal diffusivity, M b represents the Mach 
number; k*, the thermal conductivity; a 0 *, the speed of 
sound; h*, the combustion heat release; and w*, the 
reaction rate, whose dimensionless form is 


■ 

w * BzT 5 


P 

T 


n . 

Y n exp[ -E/T] 


(3.7) 


with z denoting the oxidizer-fuel ratio; n, the order of 
chemical reaction; E, the activation energy given by E * 
E*/RT 0 *; and B, the rate constant. The superscript * 
represents dimensional quantities and subscript zero gives 
the mean value at the flame edge. 


3.2.2 Solid Phase 

The solid propellant is assumed to be homogeneous, 
with no condensed phase chemical reaction. The 
dimensionless form of the heat transfer equation in the 
solid phase is 

3T 8 

ft + fiUjT^i - fT 9/ii = 0 (3.8a) 

8t 

Furthermore, assuming that the heat transfer is one- 
dimensional gives 
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3T, 3T 8 3 2 T s 

0 + r r - 0 (3.8b) 

3t 3y 3y 2 

where 

B » Ps*/Po* ' * - 3c 8 *C p */k*C 8 * 

r - r*/r* with r* - p 0 * v oVp 8 * 

Here, r denotes the burning rate at the solid surface and 
subscript s represents the solid phase. The decomposition 
process of the solid propellant at the surface is assumed 
to follow an Arrhenius law; thus. 



' 

' 1 

1 1 


r = exp 

-E. 

‘ T s 

T 

■*■3 ' 



in which E s is the dimensionless surface activation energy, 
E s =* E s */RT 0 * , and T s is the mean temperature at the 
surface. The solid-gas interface boundary conditions are 
determined by the dimensionless mass and energy balances 
across the interface, such that 


' 3T ' 

l 

' 3T ’ 

, ay 

♦ * ^ 

i 9y - 


+ rL 


(3.10) 


with £ = k*/k s * and L = (H + * - H_*)/C p *T 0 *. Here, L is the 
latent heat of vaporization of the propellant and H* 
denotes the enthalpy changes. The subscripts + and - 
represent the gas and solid side at the interface, 
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respectively. 

3.2.3 Perturbation Expansions 

When a small pressure disturbance occurs in the 
combustion chamber, every field variable will be disturbed 
from its steady-state value and can be expressed in the 
form, 

F = F <0> + eF (l> + f 2 F (2) + ... (3.11) 

where F * {p, u L , T, Y, P} and e represents the 
perturbation parameter. The superscripts in the 
parentheses indicate the perturbation order. Furthermore, 
assuming sinusoidal fluctuation of pressure with time 
renders the variables in a different form: 

F (I) = F lI, e iIut , I = 1,2,... (3.12) 

It is important to recognize that the sources in the 
second order consist of inhomogeneous terms that describe 
the nonlinearities in terms of the product of two first 
order variables. Considering the physical quantity of F ll) 
leads to the separation of each inhomogeneous term into a 
time- independent term and a term oscillating at the 
frequency of the second harmonic. Because of the fact that 
only the real parts of F^ 1 ’ and F^ 1 ’ are the physical 
quantities, the product indicates, in reality, that 
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Re(F : 1 1 * ) (F 2 
* Re[(F lr 
+ iF 2 ) (cos 


{1) ) = ReCFi 1 1 ’e 1 wt ]Re[F 2 ( 1 ’e 1 wt ] 
+ iFj) (cos wt + i sin wt)]Re[(F 2r 
wt + i sin wt) ] 


+ 


1 1 

88 - <Fi r F 2r + F X F 2 ) + " ( F l r F 2 r " F 1 F 2 )COS wt 
2 2 

1 

(F lr F 2 + F 1 F 2r )sin wt 

2 


The right-hand side of the above equation equals the real 
part of (1/2) (F : ( 1 ,+ F 2 ( 1 ’ + F 1 (1, F 2 ll> ), where the dagger 
represents the complex conjugate. Thus, the product of 
Fj* 11 and F 2 ‘ 1 1 may be replaced by 


f 1 (1, f 2 <1 > = - (F 1 ' 1, +F 2 ,1) 
2 


1 


2 


(F^’tVl’ 


+ Fj 1 1 j F 2 1 1 * ) 

+ F : ( l, F 2 e i2wt ) 


(3.13) 


and each product of two first order terms can be replaced 
by a similar expression. Since each nonlinearity yields a 
constant term and a term oscillating at the frequency of 
the second harmonic, the dependent variable F 2 ,n may be 
rewritten as 

f < 2> _ + F D <2) e i2wt (3.14) 

From the expression of Eq. (3.14), those equations for the 
second order system can be separated into two individual 
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equations, one for the time- independent equations and the 
other for the time- dependent equations. All the 
coefficients of the time- independent equation set are 
exactly the same as those of the first order set, except 
the inhomogeneous terms which consist of each product of 
two first order variables. For a higher order, the same 
argument is applicable. 

Substituting Eqs. (3 . 11) - (3 . 14) into Eqs. (3.1)-(3.5) 
and rearranging separately in the order of perturbation 
yield the final form of the governing equations 
corresponding to each order. 

Steady-State Governing Equations 

The one-dimensional steady-state governing equations 
are given as follows: 

Continuity 

p io) v io) = x (3.15) 


Energy 

3T ( 0 1 
3y 

Species 
3Y ( 0 ' 


3 2 T* 0 1 
3y 2 


a z Y (0) 

3y 2 


w (0, h 


t o ) 


(3.16) 


3y 


= -w 


(3.17) 
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State 

p (0) T (0) „ x (3.18) 

Note that the uniform pressure is retained throughout the 
flame thickness (P (0) = 1), and from Eqs. (3.6) and (3.11), 
both the dependent variables, p i0) and T (0> , are equal to 
unity at the flame edge, which is far from the origin on 
the scale of j >*.- From Eqs. (3.15) and (3.18), we have 

v 10 ’ =* T <0> (3.19) 

and from Eqs. (3.16) and (3.17), Y <0> can be expressed as 


1 

Y <0, =-(l-T ,0) ) (3. 20) 

h 


Therefore, 




2 

exp [ -E/T ( 0 ) ] 


(3.21) 


where a = 0 and the second order chemical reaction (n = 2) 
is assumed. 

Now, the equations are expressed in terms of 
temperature in the steady-state; thus, only the solution of 
Eq. (3.16) is required. Flandro [9] suggests a simple 
analytical model of the mean flow field to facilitate 
multi-dimensional analysis in the higher order system. 

This model is given in the form 
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Ui 10 ’ - u c [1 - exp(-y/R c ) ] 1 + v ,0, 3 (3.22) 

Here, u c describes the flow speed along the local 
streamline and R c is a dimensionless distance from the 
solid surface. The following boundary conditions are used 
in the steady-state solution: 

At the flame edge, 

Y (0> » o , (3.23) 

T ,0) = l (3.24a) 

or 

dT ( 01 

= o (3.24b) 

dy 

At the solid phase, Eq. (3.8) gives 

T s <0) - (T s - T c )e y/ * + T c (3.25) 

where T c is the propellant cold side temperature. The 
continuous temperature condition at the interface requires 
that 

T c o > » T < o > _ (3.26) 

8 - + 8 

The matching condition across the interface can be obtained 
by substituting Eq. (3.25) into Eq. (3.10), which yields 
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dT l0> 

dy 


T, - T c 
+ L 

n 


(3.27) 


In solving Eq. (3.16), note that a correct rate constant B 
in Eq. (3.21) has to be determined such that the system 
satisfies the boundary conditions at the flame edge as well 
as at the interface. 

Higher Order Governing Eguations 

The higher order governing equations can be expressed 
in a single form because only the source terms are 
different from each other. Presenting the source terms on 
the right-hand side of the equations in terms of G 8 , the 
governing equations are represented as follows: 
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and the reaction -rate is given by 
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Here, G s is given as follows: for the first order system 
(1=1), G = 0; for the second order system (1=2), 
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Note that, including the pressure coupling, the velocity 
coupling is significant in the source terms (Appendix B) . 
At the solid phase. Eg. (3.8) gives 
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Linearizing Eg. (3.9) results in 
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Substituting Eg. (3.36) into Eg. (3.35) and solving 
analytically yield 
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Equation (3.37) gives a boundary condition at the surface; 
other conditions are as follows: 

At the flame edge, 

Y lI> =0 (3.38) 

Assumption of an isentropic flow near the flame edge gives 
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the temperature conditions. Noting that the steady-state 
temperature gradient at the flame edge is almost zero, this 
condition can be expressed in the form 
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with 


Since the flux of each reactant species is always a 
fixed fraction of the total flux, the fuel mass flux 
fraction m f can be derived from Eq. (3.4) as 

pU i Y /i - Y #ii - 0 (3.40) 

For a one-dimensional expression, m f is given as 
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where m is the mass flux equal to pv, and assumption of the 
constant burning rate at an instant has been used. At the 
surface. 
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The normal velocity component is derived from the mass 
balance condition at the interface such that 
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Moreover, the parallel velocity component may be obtained 
using the Taylor series expansion about the surface where 
the no-slip condition must be valid. Thus, 
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Note that G 7 - G 14 are valid only for 1 = 2. For higher 
order systems, Eqs. (3 . 37) - (3 . 39) and Eqs. (3 . 43) - (3 . 45) 
are used as boundary conditions to solve Eqs. (3.28)- 
(3.32). It is necessary to have more conditions for the 
density at both sides for better solutions, and the 
pressure fluctuation has to be forced at the flame edge. 
In the case of the second order time-independent system, 
care must be taken to use boundary equations (3.35) and 
(3.39). See Appendix C for details. 
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3.3 Finite Element Applications 

The finite element method (FEM) is an approximate 
numerical procedure for solving partial differential 
equations of boundary and initial value problems. The 
basic idea resides in the application of variational 
principles or equivalent concepts such as weighted residual 
methods. It is well known that the most general approach 
to the finite element analysis is the Galerkin method, 
which is one of the methods of weighted residuals. In this 
scheme, test functions are the same as trial functions, 
known more popularly as interpolation functions. 

Typically, any scalar variable X is expressed as a linear 
combination of interpolation function and the nodal 
values of such that 

X * *> 5 X 6 (3.46) 

where B denotes the global nodes. The Galerkin finite 
element equations result from the inner product of the 
residual of a governing equation and the test function such 
as 
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Integrating this equation by parts gives rise to the 
simultaneous algebraic equations. The Gaussian quadrature 
is employed for the integrations of the spatial finite 
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element domain. Consequently, the global finite element 
equations are of the form 

[iIwA aJi + B a0 ]X B - G a (3.48) 

in which i denotes the imaginary part of the complex 
number, X fl represents the solution vector, X fi =* [p fi , u fli , 

T B' V ' and G<x ^ n ^ OItlo ^ eneous source term valid 

for the second order. 

The algebraic finite element equations are modified 
for the boundary conditions. Implementation of complicated 
boundary conditions is achieved by adopting the Lagrange 
multipliers method. The solution variables included in the 
boundary conditions have to be provided in the global form 
of the nodal values. The boundary conditions are expressed 
in the form of the boundary matrix equations, 

q-ygXjj = b-y (3.49) 

with i = 1,2, .. . ,m, with m being the number of equations. 
Thus, the modified algebraic finite element equations have 
the final form: 

iI«A a g + B a g 

In combination with the Lagrange multipliers method, 
it is also possible to impose additional Dirichlet boundary 
conditions. Implementation of Dirichlet conditions is 
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achieved by zeroing out the rows and columns of A a ^ and 
B_ a , placing 1 at the diagonal corresponding to the 
boundary nodes with all the boundary terms subtracted from 
G a . The influence of this process is propagated to all 
other equations. 

It is seen that Eq. (3.50) lends itself to a standard 
eigenvalue problem, the solution of which determines the 
existence of excited natural frequencies and modes. Using 
Eq. (3.48) with z-ero setting of G a gives the eigenvalue w 
which is the natural frequency of the given system. 

The solution of Eq. (3.50) at a given frequency is 
obtained by imposing the Dirichlet condition of the 
pressure at the flame edge. The finite element formulation 
contains two different kinds of test function used to 
represent the volume and surface integrals in the domain. 
The total number of field variables would be reduced by one 
if the density or pressure were eliminated using the 
perfect gas law. However, the stability of the matrix 
system is doubtful. 

Before calculations, the following is expected: since 
the Mach number is generally very small, the coefficient of 
the pressure term in Eq. (3.2) dominates the system unless 
the frequency considered is large enough such that other 
coefficients containing the frequency factor become 
comparable in magnitude. Therefore, in lower frequencies, 
the pressure gradient has to be negligible, thus resulting 
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in constant pressure. On the contrary, the gradient will 
become significant in higher frequencies; this results in 
pressure variance. In the latter case, severe pressure 
changes will occur if unbounded at the solid surface. Note 
that care must be exercised in expanding Eq. (3.9) due to 
the appearance of the exponential growth effect. 

For the steady-state case, as mentioned earlier, we 
calculate the eigenvalue B in which necessary initial 
conditions are satisfied. However, as will be discussed 
later, the eigenvalue is referred from the results of 
reference [4] for this study. The Galerkin finite element 
formulation of the governing equations in each system is 
described in Appendix D. 



CHAPTER IV 


DISCUSSION 

In the steady-state, the flame zone is governed by 
Eqs. (3.15) through (3.18). When an acoustic pressure wave 
at a certain frequency is imposed on the flame zone, the 
system will react according to Eqs. (3.28) through (3.32). 
The present computatipns on homogeneous solid propellants 
are performed for an adiabatic flame with a second order 
chemical reaction [Eq. (3.7)]. Instead of solving Eqs. 
(3.15) through (3.18), distributions of the field variables 
in the steady-state are calculated from the temperature 
given in Eq. (3.16). The coefficients used in this 
computation are chosen in such a way that the results may 
be compared with those in [4,9]. In view of this, the 
following parameters are utilized: z = 1, 5=0, T e =0.15, 
T s = 0.35, f = * - 1, E = 10, E, = 4, L - 0.15, 6 = 1000, 

7 = 1.2, M b = 0.003, u c = 1.0, R e = 5.0, Pr = 1.0, and h = 
1.3. Representative dimensional parameters corresponding 
to the above dimensionless values are given in Table 1, and 
these are based on a single-base propellant [12]. 

To demonstrate the validity of the theory presented 
above, a two-dimensional rectangular geometry, shown in 
Fig. 4.1, is analyzed, in which a burning surface and flame 
edge can be established as boundaries. This configuration 
was chosen to resemble a one-dimensional behavior so that 
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Table 1 Typical value and range of parameters 


Parameter 

Typical 

Value 

Range 

Physical 

Variable 

Typical 

Value 

P S 

1000 

250-1000 

-g/cm 3 

1.5 


1 ' 

- 

g/cm 3 

1.5 x 10" 3 

T c 

0.15 

- 

°K 

300 

T s 

0.35 

- 

°K 

700 

T f 

1.0 

- 

°K 

2000 

E 

10 

4-15 

cal/gmole 

40 x 10 3 

E s 

4 

2-10 

cal/gmole 

16 x 10 3 

C P 

- 

- 

cal/g°k 

0.33 

c 8 

- 

- 

cal/g°k 

0.33 

kg 

- 

- 

cal/cm°ksec 

5 x 10‘ 4 

^8 

- 

- 

cal/cm°ksec 

5 x 10 -4 

ltl 0 

1 

- 

g/cm 2 sec 

0.4 

Po 

1 

- 

atm 

9.5 

n 

2 

0.5-2 

- 

- 


0.5 

in 

CO 

• 

o 

t 

• 

o 

- 

- 

W 

- 

10 -3 -10 2 

H z 

- 

7 

1.4 

- 

- 

- 

AH 

0.15 

0.05-0.3 

cal/gmole 

2.8 X 10 3 
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the results may be compared directly with those reported in 
the 1 iterature [4,9]. 

The computations begin with finding the natural 
frequencies of the system. The resonance frequencies of 
the system can be obtained by examining the acoustic 
admittance or response function at each natural frequency. 
As mentioned earlier, the natural frequencies of the system 
are given by the homogeneous solution of the total matrix 
equation (3.50), neglecting the boundary conditions. The 
frequency range of interest in this computation is shown to 
be between w = 10" 3 and u = 10 2 , corresponding to 
approximately 80 Hz and 8 MHz, respectively. Thus, sixteen 
different frequencies in this range are chosen to evaluate 
the present study, and these are u = 10" 4 , 5 x 10" 4 , 10" 3 , 

5 X 10" 3 , 10" 2 , 5 X 10" 2 , 10" 1 , 5 X 10" 1 , 1, 5, 10, 20, 30, 
40, 50, and 100. 

Figure 4.2 shows a typical steady-state distribution 
of the field variables and the reaction rate along the 
flame thickness. In general, these data may be obtained by 
solving the energy equation (3.16) together with the 
reaction rate in Eq. (3.21). Ideally, a fourth order 
Runge-Kutta scheme may be used for this purpose. The 
result will be used as the basis of further calculations 
for combustion instability in the unsteady-state. 

Distributions of the field variables versus frequency 
for the first order system are shown in Figs. 4. 3-4. 8. 

These results are obtained by imposing an acoustic pressure 



34 


amplitude of unity at the flame edge as the Dirichlet 
condition. Conventionally, the thickness of the burning 
zone has been assumed to be negligibly small compared with 
the wavelength of the acoustic oscillation in the 
intermediate frequency range. Thus, the uniform pressure 
is approximated throughout the domain of study and is 
assumed to vary with time only. On the contrary, since the 
oscillatory pressure is regarded as a spatially 
nonhomogeneous time-dependent source in the present study, 
it is possible to investigate the response of a specific 
solid propellant at significantly high frequencies and to 
find the response even in the long flame of a double-base 
propellant. Figure 4.3 shows variations of the pressure 
distribution versus the acoustic frequencies. It is seen 
that the amplitude remains constant for u < 10. Random 
variations of these amplitudes occur for higher 
frequencies. Here, the imaginary parts representing the 
phase shift are zero. 

Figure 4.4 demonstrates distributions of various 
component fluctuations of the first order system at « * 1. 

A pressure wave with a certain amplitude striking the solid 
surface will give rise to a response in the burning rate. 
Generally, the response of the burning rate is nonlinear 
and has a complete Fourier series form that may not be 
expressed by a single frequency component. In Fig. 4.4, 
however, the pressure is uniform in the flame zone, 
although significant variations may occur at higher 



35 


frequencies . 

The temperature amplification at the surface changes 
the burning rate in Eq. (3.36), while the velocity at the 
flame edge represents the acoustic admittance for pressure 
coupling, whose magnitude and sign indicate the instability 
of the system. Since the pressure remains constant, 
implying that the acoustic wavelength is larger than the 
flame thickness in this case, the imaginary part of the 
velocity approaches a. constant slope at the edge by the 
assumption of an isentropic condition at the flame edge. 

It is interesting to note that the trend of variations of 
temperature and velocity appears to be similar, which is 
inversely proportional to the variations of species and 
density as expected. These results are comparable to those 
of T'ien [4] and Flandro [9]. 

Density distributions for various frequencies are 
given in Fig. 4.5. At w < 0.5, changes of the 
distributions versus frequency are negligible; however, 
near u = 1 a significant reduction of the amplitude occurs, 
and then for w > 1 the amplitudes increase moderately. The 
distributions in the lower frequency region are very 
similar to those in the steady-state, from which the quasi- 
steady assumption may be deduced. Note that the reverse 
peak moves in the flame zone, which implies a change of the 
reaction zone. Close to the surface, the effect of preheat 
in the solid phase increases, resulting in a higher burning 
rate. Positive amplitudes at the interface in each 
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frequency are caused by the rise in pressure and/or the 
increase of the burning rate. 

It is interesting to see that the amplitudes decrease 
while the frequency increases up to the unity and then 
increase along the frequency, from which the flame is 
expected to react sensitively at special frequencies. The 
negative amplitude is shown most notably at « ■ 1. This 
may result from the net blowing effect, and the velocity at 
that point should be positive (see Fig. 4.8). 

The normal velocity distributions are shown in Fig. 

4.8. The profile may be classified into positive and 
negative slopes, the only exception being at « * 1. The 
latter contains most of the distributions in the lower 
frequency region. Negative amplitudes appear mostly at the 
interface, relating to an increase of the density. 

Although the pressure remains constant in lower 
frequencies, the velocity varies severely. At w > 20, the 
variation is more significant since from that frequency the 
pressure varies in the flame zone. At w = 0.01, note that 
the slope is positive even though it is a lower frequency. 
This result implies instability of the system for the 
quasi-steady case. At « * 1, a positive amplitude near the 
interface results from adjusting between the density 
increase and the higher gasification rate. As indicated in 
Fig. 4.4, the slopes of the imaginary parts near the flame 
edge are constant. 

Rearranging the real part of the velocity at the flame 
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edge gives the distribution of the acoustic admittance, 
whose magnitude and sign indicate the amplifications or 
damping ability of the flame subject to the acoustic 
disturbance. Figure 4.9 shows the curve-smoothed acoustic 
admittance and burning rate from Eq. (3.36) at each 
frequency. 

At « = 1, there is a definite increase in the burning 
rate. The burning rate increases by increasing the heat 
transfer effect due to closer movement of the reaction zone 
to the surface. There are actually several secondary 
effects such as structural change in the solid phase and 
change of the chemical kinetics in the flame. These 
effects cause the feedback to the solid phase, resulting in 
a change of the burning rate. These effects, however, are 
not considered. 

Figure 9 also reveals a resonance in the condensed 
phase near u = 0.01, indicating that the system is 
unstable. This verifies the result of Denison and Baum [2] 
and previous laboratory measurements which have shown that 
the response consists generally of a single peak in the 
range of frequency for which the quasi-steady approximation 
appears to be accurate. Some negative peaks exist at the 
other frequencies, implying that the resonance in the gas 
phase tends to damp the oscillatory motion. At most higher 
frequencies, the system seems to be unstable. 

The real part of the burning rate shows a trend 
similar to the acoustic admittance at the quasi-steady 
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region, although the magnitude is significantly different. 
But these trends differ from each other at the higher 
region. Over w * 100, the profile tends to have a second 
mode oscillation as the pressure varies in the flame zone. 

It must be emphasized that the admittance function alone is 
not sufficient to describe the stability of the system 
unless the velocity coupling is concerned. 

Temperature distributions are given in Fig. 4.6. 

Since the temperature is related to the density by Eq. 

(3.5), it can be analyzed easily by comparing the results 
with those in Fig. 4.5. Accordingly, the distribution of 
the temperature shows the profile opposed to that of the 
density. Temperature rises are predominant along the 
reaction region for w > 1. Note that no significant 
temperature changes occur at the interface as imposed by 
the boundary conditions. Negative amplitude near the flame 
edge may be caused by the stagnation phenomena in which the 
density increases. Because of the pressure unity, the 
imaginary parts of the temperature are given as reciprocals 
of those of the density. 

The fuel consumption at a given frequency is shown in 
Fig. 4.7. Trends for species distributions are opposite to 
those for temperature distributions in that negative 
maximum occurs along the reaction region. Note that 
linearly diminishing the fuel near the flame edge affects 
the temperature changes slowly (Fig. 4.6). 

The results of the second order perturbation are shown 
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in Figs. 4.10-4.28. As previously indicated, each variable 
of the second order response to acoustic disturbance has 
two components: one time-dependent component that 
oscillates at twice the fundamental frequency and one that 
is time-independent. The latter represents a shift in the 
mean value, thus causing a shift of the mean burning rate. 
It should be noted that the variations for the second order 
system may be characterized by the sum of those two parts, 
with the factor e 1 2 ut . ranging between -1 and 1 and zero 
indicating the case of time-independence. 

Figures 4 . 10-4 . 15 show the distributions in the second 
order time-independent case. The calculations are also 
conducted by imposing the pressure of unity at the flame 
edge as the Dirichlet condition. General trends show that 
shifts in the mean values are evident. 

In Fig. 4.10, the pressure varies from w = 5, which is 
half of that in the first order system. For the higher 
frequency region, the oscillatory motion is significant, 
which may affect the velocity distribution. Figure 4.11 
shows all of the variables for w = 1, indicating that 
shifts in mean values may have the affect of damping; 
however, the opposite may be true for higher frequencies as 
demonstrated in Figs. 4.12-4.15. The trends are roughly 
similar to those of the first order illustrated in Fig. 
4.4a, but they have different amplitudes because of the 
nonlinearities in the higher order system. In Figs. 4.12- 
4.15, the distributions of field variables in lower 
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frequencies appear to vary negligibly along the flame 
thickness, while the opposite is true for u > 1. In 
particular, since the pressure varies at « > 5, the 
velocity changes significantly at those frequencies, 
implying that the system is unstable at higher frequencies. 
As mentioned earlier, a shift of the mean burning rate is 
the most significant effect in the time- independent system. 
This result will be discussed in Fig. 4.28. 

The distributions in the second order time-dependent 
case are shown in Figs. 4.16-4.21. Here, the pressure also 
varies from u = 5 and oscillates for higher frequencies. 
Note that the amplitude variations are smaller than those 
in the time-independent case. The trends seem to reduce 
the total effect of the second order system. 

Figures 4.22-4.27 show the variations of field 
variables in the second order system by summing time- 
independent and time-dependent effects, with the factor 
e i2ut » i for time-dependent effects. Note that amplitude 
variations of each variable in the flame zone are 
significant at « = 1, which is the same result obtained for 
the case of the first order system. Finally, Fig. 4.28 
reveals that the shift in burning rate versus frequency is 
consistently negative, as asserted by other investigators. 
The offset is relatively small in the quasi-steady region, 
but increases with oscillatory motion along the 
frequencies. The variations of the acoustic admittance are 
similar to those of the first order system, except that 
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oscillatory deviations are more predominant at higher 
frequencies in both time-independent and time-dependent 
cases. 

Parameter studies are conducted for the first order 
and are summarized as follows. Decreasing the density 
ratio 6 affects the variables shifted slightly to the 
negative direction, keeping the distribution profiles 
constant. Changing the latent heat of solid L exerts a 
negligible effect on the variables, but a very small value 
of L shifts the system toward instability. Increasing the 
surface activation energy E 9 or the gas phase activation 
energy E reduces the magnitude of the variables, keeping 
the same profiles, changes of the rate constant and 
viscosity effect coefficient strongly affect the system, 
such that every aspect discussed herein will change. 

The present study could be extended to the multi- 
dimensional case by introducing the appropriate axial mean 
flow field. It is well recognized that fluctuation of the 
gas velocity parallel to the propellant surface affects the 
burning rate in terms of velocity coupling; therefore, this 
quantity must be considered together with pressure coupling 
for a satisfactory measure of stability. 

A simple calculation has been accomplished using the 
artificial axial flow velocity in Eq. (3.22). However, 
difficulties of the boundary conditions could not be 
eliminated. A test run shows that the existence of a small 
amount of the axial flow reduces differences between the 


variable amplitudes in each frequency considered, thus 
leading the system toward stability. 
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CHAPTER V 
CONCLUSION 

A multi-dimensional numerical model for the premixed 
flame acoustic instability is proposed and solved using the 
finite element method. The governing equations are 
perturbed to the second order and formulated with the 
Galerkin finite elements. The results have direct bearing 
on the validity of published theories of solid propellant 
combustion instability at the lower frequency region where 
the uniform pressure assumption is valid. Extended studies 
are made on the higher frequency region and some new 
results are obtained. 

Under the restricted boundary conditions, the 
following conclusions, based on numerical computations, are 
reached : 

(1) The stability characteristics for the low frequency 
range in the first order system have been verified to 
be the same as those reported in the literature. For 
example, the acoustic admittance is controlled by the 
burning rate, negative for low frequencies, whereas 
the opposite is true for high frequencies. 

(2) Instabilities are likely to occur at high frequencies 
( w > 10) . 

For the second order system, the trend is similar to 
the first order system for low frequencies, but 


( 3 ) 
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instabilities may appear at frequencies lower than 
those of the first order system. 

(4) The most significant effect of the second order system 
is that the admittance is extremely oscillatory 
between u - 1 and w - 10. 
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Fig. 3.1 Domain of study for one-dimensional case 
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Fig. 4.1 Finite element geometry. 
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Fig. 4.2 Steady-state distributions of field 
variables in flame zone. Parameters 
used in the calculations are given 
in Table 1, including Pr=l and Mb= 
0.003. Reference values for non-dim- 
ensional ization are chosen from the 
flame edge. ( p : density ,T : temperature , 
Y:species of fuel, v:velocity, P:pre- 
ssure, and w:reaction rate). 
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Fig. 4.4a First order distributions of field 
variables at cj=l(real parts). 
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Fig. 4.8 First order velocity distributions 
vs. frequency. 
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L7a Second order time-dependent distr: 
butions of field variables at to = 
(real parts) . 
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Fig. 4.18 Second order time-dependent density 
distributions vs. frequency. 
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Fig. 
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4.20 Second order time-dependent species 
r (fuel) distributions vs. frequency. 
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Fig. 4.21 Second order time-dependent velocity 
distributions vs. frequency. 
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Fig. 4.22 Second order pressure distributions 
vs. frequency. Both time-dependent 
and time - independent coefficients 
are added for denoting the maximum 
deviation from mean value. 
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4.23 Second order distributions of field 
variables at to = 1. 
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Fig. 4.28 Offset burning rate and acoustic 

admittances of the second order time 
-dependent and time-independent sys- 
tems. Negative value of the burning 
rate at the frequency region verifies 
the past investigations. 
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APPENDIX A 

DERIVATION OF NON-DIMENSIONAL GOVERNING EQUATIONS 


1. Continuity Equation 

The dimensional continuity equation is of the form 


dp* 


at" 


+ (P*u i*) ,1=0 


Substituting Eq. (3.6) yields 


Po*v 0 * d P 1 


— + — (p 0 * v o*P u i) , i = 0 

i* at a* 


Dividing (A. 2) by p 0 *v 0 */F yields 


dp 

— + (PU i) ,1=0 

at 


which is equivalent to Eq. (3.1). 

2 . Momentum Equation 

Neglecting the body forces, the dimensional governing 
equation is given by 


3uj 


at" 


+ p*n i#j u i * + P yi - u. 


u 


i / j ] 


1 * 
+ - u 
3 


J / i ) 


= 0 


Substituting Eq. (3.6) gives 


* 2 


3u 


*2 


P 0 — p + Po P u i , j u j 

j l* at jt* 



(A.l) 


(A. 2) 


(A. 3) 


(A. 4) 
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I v * 
V 0 


- u 


Ui 


Uj 


*2 1 t 1 J - 0 * 2 • 1 j 


n 


= 0 


(A. 5) 


Dividing the equation by p„*v 0 * 2 /j!* results in 


at 


+ P n ifi n i + 


u 


Po*V 0 * 2 


Pi - 


p 0 *v 0 *r 


u 


i / j j 


+ " U j,ij 


= 0 


Since 


(A. 6) 


Mb 2 = 


( \T * \ 2 


V a o J 


Po V 


*„ * 2 


7P r 


where M b is the Mach number at the reference point, then 


„ * 2 
P 0 v o 


7Mi 


(A. 7) 


For the last term, 


l* = 


P o *C p *V o * 


with a* being the thermal diffusivity. Thus, 


u 


p 0 *v 0 *r 


U*p o*C p *V 0 * 
p 0 *V 0 *k* 


U*C p * 


= Pr 


(A. 8) 


Substituting Eqs. (A. 7) and (A. 8) yields 


3Ui 1 

p + 7 p ,i - Pr 


at 


7Ml 


U i / j i + _ U j / i j 


1 

3 


= 0 


(A. 9) 
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which is equivalent to Eq. (3.2). 


3. Enercrv Equation 

Neglecting the radiation effect, the dimensional governing 
equation is given by 


3T’ 


3P’ 


n 


P*c p* — + P*C p*T #i *Ui* 


k T j i - l w a h a = 0 (A. 10) 

3 1 a=l 


Substituting Eq. (3.6) yields 


1 3-T v *T * 

p 0 *C p *V 0 *T 0 * — p — + P 0 *C p * — — ■ 

St* at jt* 


P*, iUj - 


P 0 *v 0 * 3P 


r at 


- k* T i , - 

0*2 ' 


* 2 

Po v o 


C p *T 0 * l W a h tt = 0 


(A. 11) 


Dividing Eq. (A. 11) by p 0 C p v 0 *T 0 /% gives 


3T 

P — + P T , i u i “ 


at ' p 0 *T 0 *c p * at 


3P 

T, i i - S w a h a = 0 

a 


(A. 12) 


Since P 0 * - P 0 *RT 0 * , 


P 0 * R C v * 7-1 

Po*T 0 *C p * C p * C p * 7 

where R * C p * - C v * and 7 = C p */C v * are used. Therefore, the 
non-dimensional governing equation is of the form 


3T 7 - 1 3P 

P — + pT jUj T h ~ Wh = 0 (A. 13) 

3t 7 3t 


in which the last term is obtained from Eq. (A. 30). Equation 
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(A. 13) is equivalent to Eq. (3.3). 

4. Species Equation 

The dimensional governing equation is of the form 


ay. 


at' 


+ p Y a) jUj — y a # i i ** W a — 0 , a 1/2 / ••• / 


where N denotes the total number of species considered. 
Substituting Eq. (3.6) yields 


v 0 * 3Y a v 0 * k* 1 

Po „ * P _ + P 0 . * P Y a,i u i " _ * 2 Y a,ii 


r at 


r 

* V * 2 
P o v o 


c D * J 2 


w a - o 


a 


Dividing Eq. (A. 15) by /s 0 *v 0 */^* gives the following 


ay. 


+ p y tt/i Ui - 


at c p *j!*Po*Vo* 


Y a/ii - 


v 0 *r 


w a = 0 


Since 


*y *P * 

Po V 0 V 'P 


The coefficient of the third and the last terms are to be 
Therefore, 


ay. 


at 


+ P Y a,i u i - Y a,ii - y a w a = °> a * f *°/P 


N 

(A. 14) 


(A. 15) 


(A. 16) 


unity. 


(A. 17) 
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where y a is defined by 


W a ('a" - v a . ' ) 


W ft' 


(A. 18) 


The three mass fraction equations for the fuel, oxidizer, and 
product in Eq. (A. 17) are similar to one another, allowing 
treatment of the fuel equation alone with the following relation 


— - — » constant • 
v a v f 


(A. 19) 


Therefore, Eq. (A. 17) for the fuel only is enough to describe the 
whole system and is given as 


3 Y 

p + pY f i u i - Y f/ii + w f = 0 (A. 20) 

at 


This is equivalent to Eq. (3.4). 

5. Reaction Rate of Fuel Species 

For a one-step forward chemical reaction, 


K f '[f] + K 0 »[0] ■> V'[P] (A. 21) 

where v denotes the stoichiometric coefficient, and f,o,p 
represent fuel, oxidizer, and product, respectively. Using the 
Arrhenius law, the reaction rate of the fuel can be expressed by 


w 




r p* 

n N 

TT 

' h 1 

, RT* , 

j=i 

, Wj , 


e -< E * / R T * ) 


(A. 22) 


in which B* is the frequency factor for the gas phase; R, the gas 
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constant; Wj , the molecular weight of the species; E*, the 
activiation energy; and 5 and n are constants. For the oxidizer 
and products, the reaction rate can be written as follows: 

y It - y » 

w a * - — — w* , a = o,p (A. 23) 

' 

The the last term of Eq. (A. 12) becomes 

- I w a *h a ° = w*h* (A. 24) 

where 

* 2 W«(r«' “ ^')h a 0 

h* = 

W f K f ' 

is the heat of combustion per unit mass of fuel consumed. Using 
the definition in Eq. (A. 18), 


".('a" - y a ' ) 

v a " 

W f v t ' 

and the following relation in Eq. (A. 19) 


(A. 25) 


Ya Yf 

— - — - constant 
•'a v f 


(A. 26) 


the reaction rate term can be non-dimensionalized. 

Since both the fuel and oxidizer must vanish at the flame 

edge, 


Y 0 Yf 

=0 (A. 27) 
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Therefore, 


Y 0 = — Y f (A. 28) 

v ( 

From Eq. (A. 25), v f = -1; thus, 

y 0 =* -y 0 *f 

or 


Wo'o’ 

Y o - Y f (A. 29) 

W f K f 


The non-dimensional reaction rate can then be written as 


w 


BzT s 

i 


P 

T 


Yp n e 


E/T 


(A. 30) 


where 


W 0 'o' 

2 ~ 

B*k*T 0 * J p 0 * n 
B . 

Po* 2 v 0 * 2 c p n Wj n j 

j=i 

E* 

E = 

RT 0 * 

and this is equivalent to Eq. (3.7). 
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Heat Transfer Equation in Solid Phase 

The dimensional governing equation is given as 


* _ * 
Pa c s 


+ Ps * c s *T* /i u i * - k s *T* #ii = 0 


(A. 31) 


with c s * and k s * being the specific heat and the thermal 
conductivity of the propellant, respectively, and y = 0 is 
attached to the surface. Substituting Eq. (3.6) yields 

v 0 * 9T, . T 0 * T„ * 2 

P ,*c s *T 0 * — > P 3 *c s * — v 0 *T s<i Ui - k s * — - T s j j = 0 

r at i* a * 2 

(A. 32) 

Dividing Eq. (A. 32) by p 0 *c s *T 0 *v 0 */j>* results in 


3 + J3T s t j u j — 


V *rp * 
J '-S 1 0 


a j. * o * 

p 0 c s v 0 l 


T S/ii - 0 


(A. 33) 


Since j?* = k*//s 0 *v 0 *C p * , the last term can be written as 


v *<v * 

^ S x 0 


v */-■ * 
'■'p 


P 0*Cg*V o *r k*Cg* 


(A. 34) 


Therefore, Eq. (A. 33) becomes 


6 + BT 8 , jUj - (T s j j = 0 

3t 


(A. 35) 


In the one-dimensional case, since p 0 *v 0 * = p s *r* if r is defined 
by r = r*/r* (actually r* is the same as Uj* in Eq. (A. 31)), then 


(A. 36) 
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which is equivalent to Eq. (3.8). 


7. Solid-Gas Interface Boundary Condition 

The surface energy balance across the interface can be 
expressed as follows: 


-k 


9T’ 


9y* 


+ p*V*H + * = -k, 


f 9T* ) 
9y* 


+ P *r*H_ 


(A. 37) 


with subscripts + and - bjeing used to indicate the gas and solid 
side of the interface, respectively, and H* denotes the enthalpy. 
Using Eq. (3.6) and defining £ = k/k s * and L = (H + * - H_*)/C p *T 0 * 
yield 


T 0 * 
k* 

' 9T 1 

+ p s *r*C p *T 0 *rL = 

m * 
* 0 
-k s * * 

' 9T ' 

r 

i 9y > 

+ 


, 9y , 


or 


k* f 3T ^ Ps *r*c p *j>* 


*s* 


ay J 


rL = - 


9T 


9y 


(A. 38) 


Since p s *r* = p 0 *v 0 * and j>* = k*/p 0 *v 0 *C p * , then 




' 9T 
9y ) 


+ £rL 


9T 

9y 


or 


r 9T 

l 

' 9T ’ 

. ay , 

♦ « 

i 9y , 


+ rL 


which is equivalent to Eq. (3.10). 


(A. 39) 



APPENDIX B 


PERTURBATION OF MODEL EQUATIONS 

1. Perturbation of Reaction Rate 

The dimensionless reaction rate equation is given by Eq. 
(3.7). After taking 5*0 and n = 2, we have 


,< o ) 


w ,v " + €W (1> + f 2 w (2> = Bz[P 10 ' + eP 


( o > 


,( i ) 


2 t i 2 rml 0 I 


+ + € T <n + £ 2 T ,2, ] -2 [Y (0> + eY 

j. , 2 v ( 2 ) -i , 


v(l) 


■ - e ' 

/ t ( i > 

*1 1 - i 

ip ( 2 ) 

f 2 

- l • 

T ( 0 ) 

T {0> 

T ( 0 » J 



t 0 ) x 2 


Bz 


,( o ) 


y ( 0 ) 2 e -E/T 


{ 0 ) 


2P 


( 1 ) 


[ 2P 


( 2 ) 


1 + f 


,( 0 ) 


+ £ 


[ P 


( 0 ) 


.( 1)2 


.( 0)2 


+ 0(f 3 ) 


2T 


( l ) 


1 - £ 


, ( 0 ) 


- £ 2 


2 T 


( 2 ) 


3T 


( 1)2 


,( 0 ) 


.( 0)2 


+ 0 ( € 3 ) 


+ 0 ( £ 3 ) 



r 2Y <1) 

1 + r 1 

-2 

f 2Y (2> 

Y< 1 > 2 

_ . 


JL i C T 

Y (0) 

C 

Y ,0> 

Y 10)2 


exp 

E 

r T <1> 

c 2 

f T ( 2 ’ 

.j. ( 1 ) 2 

T < 0 ) 

T ,0) 

t 

T ,0 > 

T ( 0)2 


+ 0 ( £ ’ 


= W <0, [1 - 2£T t ' - £ 2 (2T 2 ' - 3T X ' 2 ) + 2 £P X • + £ 2 (2P; 


+ P^ 2 ) - 4£ 2 T 1 'P 1 '][1 + 2 £ Y j 1 + £ 2 (2Y 2 ' + Y x ' 2 ) ] 


( E 


+ £ 


. ( 0 ) 


T, ' + - £ 


\T ( 0) / 


T i 2 


1 + £' 


.( 0 > 


(T 2 * - T j 


■ w (0) [l + €(-21!* + 2P! 1 ) + £ 2 (-2T 2 ’ + 3T! 12 + 2P 
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+ Pj' 2 - 4T X 'P X ')][1 + 2 fY x ' + e 2 (2Y 2 ' + Y x ' 2 )] 


E 1 [ E ' 2 


+ e 


■ ( o ) 


T, • + - r 


,( o > 


T x ' 2 + e 2 


i ( o ) 


(t 2 * - iv 2 ) 


w 


( 0 ) 


1 + e(-2T x ' + 2P X ') + e 2 (-2T 2 ' + 3T X * 2 + 2P. 



E 


w 

+ 

H 

O4 

1 

2Y X 1 + T x ' 

T (0> 1 J 

+ < 2 


2Y X ' 


I ( 0 > 


(-2T X ' + 2P X ' ) + 2 f ' 


— Y x 'T x 1 + € z (2Y 2 1 


+ Y x • 2 ) + - e 2 
2 


E ) 


i ( o ) 


+ ( 


.( 0 ) 


(T 2 * - IV 2 ) 


w 


( 0 ) 


1 + € 


-2T X ' + 2P X * + 2Y X ' + 


,( o ) 


+ € 2 (“2T 2 ' + 3T x ' 2 + 2P 2 ' + P x ' 2 - 4T X 'P X ' - 4Y l 'T l ' 


+ 4Y X 'P x • - 


2E 


i ( 0 ) 


i 2 


2E 


it 0 ) 


T x 'P x ' + 


2E 


, ( o ) 


Y x ' T x ' + 2Y 2 • 


+ Y x ' 2 + - 
2 


( E 


l T 


t o > 


T, ' 2 + 


E E 

T,' - 


1 1 o ) 


it o > 


T x ' 2 


(B.l) 


where F x • = * i > /F * o » and Fj , , ^.t 2 » /p t o > are used> Therefor6/ 

( P l0) \ 2 


W l0 » = BZ 


A I 

w' 


w< 2> - w ,0 » 


it 0 ) 


yt 0 ) 2g-E/T 


t 0 ) 


= w « o ) 


r E 

ipt l ) 

2 yt 1 ) 

2 p (1) - 

. 

T t0> J 

T t 0 ) 

y ( 0 ) 

1 

P« 0 ) 


it 0 ) 


\ t‘ 2 * 


- 2 


2Y 


( 2 ) 


it 0 > y t 0 > 


2P‘ 2 ’ 


,t o ) 


(B.2) 


(B. 3 ) 


+ 3 


1 - 


it 0 ) 
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6T 


( 0)2 


' p ( 1 ) 

2 

y ( 1 > 2 

£>< 1 > 2 

' 2E 

- 4 

p< 1 ) 

p< 1 > 

, T* 77 

2 

y ( 0 ) 2 

+ p< 0 ) 2 4 

T (0 > 

p< 0 > 

p( 0 ) 


r 2 E 

' 

- 4 

y 

p( i ) 

y ( 1 ) 

p< 1 ) 

y( 1 ) - 

T ( 0 ) 

p ( 0 ) 

y ( 0 > 

T 

p( o ) 

y ( 0 ) 


(B. 4 ) 


Since 


P 

and 


( i > 


p (0) T (l) + p \ i » T 


< 1 >m< 0 > 


£< 1 > 


,( 0 ) 


,< 1 ) T ( 0 ) 


(B. 5 ) 


,( 2 ) 


p ( 0>T ( 2 > + p 1 1 ' T ' 1 ' + p' * 'T 


( l >m< 1 ) 


< 2 >m( 0 ) 


£< 2 ) 
m < 0 ) 


+ P 


( 1 > £,< 1 > 


+ P 


( 2 >m< 0 ) 


we have another form of the reaction rate such that 


w 


( l ) = w ( o > 


w< 2 > = w <0) 


( o ) 


( 0 ) 


( 1 ) 


p (2 > + 


,< 0)2 


,( 0)2 


,( 1 > 


T< 2 > + 


•10) 


• ( 0 ) 


Y ( l ) 


y‘ 2 * + 


( 1)2 


< 0)2 


E 

r E 

- 1 

J 

p ( 1 > 2 


y( 1 ) 

2 2E 

A ( 1 ) 
P 

1 ) 

p( 0 ) 

, 2T <0> 

p ( 0 > 2 

r 

Y 1 0 1 

2 p< 0 ) 

P <0 ’ 

p< 0 ) 

A ( 1 ) <r( 1 ) 

P x 

2E T ( 1 1 

A f 

y‘ 

“ 1 





+ 4 


( 0 ) v< 0 > 
p * 


ip( 0 ) ip( 0 ) y ( 0 ) 


(B.6) 


(B.7) 


If we substitute T (0> for l/p ((n and (1/h) (1 - T (0> ) for Y 101 , 
then 


w ll > = w <0) 


2T 


(0) A I1) 


,( 0)2 


T< 1 } + 


2h 


1 - T 


( o ) 


y < 1 > 


(B.8) 
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w ,2) - w‘°> 


2T <0, p ,2> + 


,( 0)2 


,( 2 ) 


2h 


1 - T (0> 


Y< 2 > + ip( 0 ) 2* ( 1 ) 2 


,p( 0 ) 3 
2E 

m( 0 ) 


2T 


( 0 ) 


- 1 


,( 1)2 


(1 - T l 0 * ) 2 


A A 4hT t0) 

p ( 1 > T ( 1 > + — + 


Y< l > 2 


2hE 


® ^ 1 — T ^ 3 ^)T^ 3 ^ 2 

Combining Eqs. (B.6) and (B.7) yields Eq. (3.33). 


.ip C 1 I y 

( B . 9 ) 


2 . Perturbation of Burning Rate 

From Eq. (3.9), the burning rate r can be perturbed in the 

form 


.( o > 


r"" + cr ( 1 ’ + c^r 1 *' = exp[-E s (T 8 

+ e 2 T 8 , 2 , )' 1 ]e E s/T s = e E 8 /T 8 exp) 


2 *( 2 ) _ 


( 0 ) 


+ eT 


-E, 


( o ) 


( l ) 


( l ) 


1 - £ 


( 0 ) 


- e 2 


( T 


( 2 ) 


( 1)2 


l T s 1 0 ’ T s t0)2 


e E s /T s e“ E s ' T s ° * expl 


( i ) 


1 s s 

( 


«p g ( 0 ) 1J, ( 0 ) 


i T. 


( 2 ) 


( 1)2 


( 0)2 


= e E s /T s e" E s /T s 0) 


exp| 

E 


( o ) 


( o ) 


1 + £ 


s T s 


( l ) 


m ( 0 ) m ( 0 ) 
*8 A S 


+ £ ' 


E s 

( m (2) 
A S 

T s ‘ 1 » 2> 

1 

1 

f E s 

m ( 1 ) 
A S 


+ 0(£ 3 ) 

m ( 0 ) 
A S 

m ( 0 ) 
V X S 

m (0)2 

T s ) 

T — 
2 

m < 0 > 
v A S 

T ( 0 ) 
s * 

1 


(B. 10) 


Therefore, 


r ( 0 ’ - exp 


— “P 

l 

1 



m ( 0 ) 
1 S 

T s 



= l 


(B. 11) 
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where T # <05 *T # at y = 0 is used, 


t 1 > „ _( 0 ) 


( 0)2 


T s (1) = r <0, c 1 T s ,l) 


£( 2 > 


,( 0 ) 


( 0)2 


' A T s 

T s ‘ 2> 


( 1)2 


' S A 


( 1)2 


( 0 ) 


r t0, Cl 


T g 1 2 ) + 


( c. 


■ s 

1 ' 


2 T 


2T, 


( 1)2 


( 0)2 8 


Combining Eqs. (B.12). and* (B.13) yields Eq. (3.36). 


(B. 12) 


(B.13) 
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APPENDIX C 

DERIVATIONS OF BOUNDARY CONDITIONS 

1. Temperature Condition at Interface 
(a) First Order Condition 

The first order solid phase energy equation (3.35) is 
rewritten as 

dr A 3T S ( 0 ’ 3 2 r 

iwflr + r <0) — + r‘ 1> C - 0 (C.l) 

ay ^ ay ay 2 

where r denotes T s (1> for convenience. Knowing that r (0> = 1 and 
aT 8 ,0, /3y = (1/f) (T s - T c )e y/ ^ from Eq. (3.25), we have 

9 2 t d-r A 

f ilwBT = c 2 r ll) e y/f (0.2) 

3y 2 dy 

in which 

1 

c 2 = - (T s - T c ) 

C 

The general solution of Eq. (C.2) has the form 

r = r h + T p (C. 3 ) 

with subscripts h and p denoting homogeneous and particular 
solution, respectively. The characteristic equation is then 
given by 

n 2 - \ - iuB = 0 (C. 4) 
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Thus, 


! - — (1 + -J l + 4i» 
2 ( 


) 


X x ' = — (1 - -J 1 - 4i«flf ) 
2 { 


(C.5) 


which gives the homogeneous solution such as 


T h “ A 0 e X lY 


(C.6) 


where r h ■» 0 at deep-in solid (y -» -®) is considered. 
The particular solution T p is of the form 


r = A 1 eY /i ‘ 


(C.7) 


Since dT p /dy - (l/f)r p and d 2 r /dy 2 » (l/f 2 )r p , we have 


*! - - 


( 1 ) 

r* 1 'c. 


iufl 


(C.8) 


Therefore, from Eqs. (C.3), (C.6), (C.7), and (C.8), r is given 
by 


r ( 1 ’c 

r = A 0 e x i y e Y/ ^ 

iwfl 


(C.9) 


From the fact that r = r s at y = 0, 



(C.10) 


which yields 
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,( l ) 


T 8 + 


C 2 ' 


iuB 


.ill, 


x i y - 


.y/i 1 


iufl 


(c 


and 


9 r 

3y 


T S + 


r 1 1 > c 2 ' 


i uB 


.m, 


X lY - 


i uBf 


e 


Y/f 


(C 


At y * 0, 


ar 

ay 


y a o 


T s + 


r“’o, ' 


iwB 


.( i > 


iwflf 


(c 


Substituting Eq. (C.13) into Eq. (3.10) yields 


r aT 


( 1 ) 


1 

r a t ’ 

♦ i 

. 3y , 


y = o 


( l > 


T| ,1> + 


+ r (1, L 


r (1, c, n 


* ( 1 1 
r 1 1 ' c. 


iwfl 

C 1 C 2 A 

— T, 1 " 
iwB 


iwBf £ 


+ r ( 1 ’ L 


c,c, 

^ * A / < t A 

T a 1 ' + C X LT S 


iufif £ 


(C 


where r‘ 1 » * r ( 0 » (E 8 /T c 2 ) T 8 ( 1 » « c^^ 1 ' is used. Therefore, 


3T 


( l > 


T + (1 » = 


3y 


x l 

\ c lC 2 

1 4- 

CjCj 

— 1 — T 

. % 

X T 

iwB , 

• C j Jj 

iwflf £ 


(C 


. 11 ) 


. 12 ) 


.13) 


( 1 ) 
.14) 


.15) 
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(b) Second Order Condition 

The second order solid phase energy equation is rewritten as 


i2w5r + r‘ 0 1 — + r ( 2 ’ 
ay 


3T, 


< o ) 


3 2 r 


ay 


- r = -r (1> 

ay 2 


3T, 


( l ) 


ay 


or 


a 2 r 
f — 
3y 2 


3t 

ay 


i2uflr * r 


( i ) 


A 

3T , 


( l ) 


+ r 


( 2 ) 


3T 


( o ) 


ay 


3y 


(C. 16) 


With the same argument of the first order solution, the 
characteristic equation is given by 

n 2 - X - i2wS - 0 (C. 17) 


of which solution is 


X 2 
X 2 


= — (1 + -J l + siwfif ) 
2 £ 

! I 


I = 


(1 - 


1 - 8iuB£ ) 


2 f 


(C. 18 ) 


Therefore, the homogeneous solution is 

r h =* B 0 e x 2 Y (C. 19) 

Substituting Eq. (C.ll) into Eq. (C.16) yields the particular 
solution to be of the form 

T p * B 1 e x i y + B 2 e y/f (C. 20) 


Since 
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3 T n B 2 

- \ 1 B 1 e x i y + — e y/ ^ 

3y f 


and 


3 z r, 


ay 


B- 


X 1 2 B 1 e x i y + — e y/f 


we have, from Eq. (C.16), the following equations: 


For the coefficient of the e x i y term, 


+ XjB! + F - i'X 1 2 B 1 = 0 


(C. 21) 


where 


F - r (l, X 


T e 1 1 * + 


r 1 1 1 c, 1 


iwfl / 


For the coefficient of the e y/ ^ term, 


B 2 A S 

i2 uBB 2 + — + r tz, c 2 - — B, - 
f f 2 


r { 1 * 2 c ; 
iwfii - 


= o 


(C. 22) 


Thus, 


B, = 


B, = 


^Xj 2 “ X ^ — i2 uB 


-r (2> c 2 + r 


( 1)2 




i2wB 


(C. 23) 


and 
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T ■ T h + T P ■ B 0 e* 2 Y + 


r x x 2 - X 2 - i2u3 


X lY 


-r <2> c 2 + r <1)2 


Cj/iuBi 1 


i2 wfl 


(C . 24 ) 


Since r = r s at y = 0, then 


B, 


-r (2, c 2 + r ( 1 ) 2 c 2 /iwflf 


T s " 


n x 2 - X 2 - i2ufi 


i2 u& 


(C. 25) 


which yields 


-r ( 2 , c 2 + r l,, 2 c,/i«« ) 


T s - 


rx 2 2 - x 2 - i2ufl 
F 


i2uB 

A 


, X 2 y 


-r (2, c 2 + r‘ 1 ’ 2 c 2 /iw6f 

e^i y + e y/ ^ 

fXj 2 - x 2 - i 2 ufl i 2 uR (C . 26 ) 


and 


3t 

ay 


-r <2, c 2 + r l 1 ' z c 2 /iu6r \ 

— ■ a x 


( 1 ) 2 , 


U “ 


f X 2 2 — X 2 — i2 wJ3 


X t F 


-r l 2 *c 2 + r ‘ 1 1 z c 2 /iw Rt 

e x i y + e y/f 

fX x 2 - Xj - i2wfl i2 uRi (C. 27) 


i2w3 
.(1)2, 


e x 2 y 


Substituting Eq. (C.27) into Eq. (3.10) yields 


3T ,2) ] 

i 

1 H 

r a r ' 

. ay > 


l ay > 


+ r ,2) L 


y = 0 


T. 1 21 - 


r ( 2 ] c 2 - r l 1 ’ 2 c 2 /iuBf 


n t 2 - X x - i2wJ3 


i2 ufl 
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1 

+ - 


XiF 


( 2 ) r _ £( 1 > 2 


r ' ‘ ' c 


l n x 2 - X x - i2uB 


2 - r ' 11 Cj/iwBf 
i2wfirj 


Since 


+ r (2, L 

(C. 28) 


r ‘ 2 ’ = 9 


1 

f 1 E, 

t, <2> + — 

I * 

L T s 

i 2 T s 


( 1)2 


T s (2> + 


f C, IX 


( 1)2 


t 9 T l 2 * ] 

L 

| X 2 

A . 

T . 

' 2 T, / J 

F c c T ( 2 ) 

( 2 ) - * , ° 1 C 2 T S 

. ay ; 

■U/ 

1 

A 5 ' 

£ X x 2 - X x — i 2 uS i 2 ufl 


f C 1 

1 1 

I 


C 1 C 2 


T s ) 


S ' A 


.( 1 ) 2 , 


l 2 uB 


s 


C C T < 2 > 
C 1 C 2 1 8 


CiC 2 


( 1)2 + 

Cl 1 
2 


2 w 2 B 2 r 


i 

+ - 


XiF 


t r x ! 2 - x x - i2us 


T 8 ) A 


.( 1 ) 2 , 


i2wfif ^ 


( 1)2 


i 2 uBn 


2 w 2 fl 2 f 2 £ 


+ c x 

T 8 ‘ 2 » + 

'Cl 1 ’ 

T s ( 1 > 2 


. 

.2 T, , 



(C. 29 ) 


Therefore, 


( 2 ) 


' 

'aT <2> 


3 y - 


X,F 


X 2 c 1 c 2 


fcj 1 X 


TJ 


+ - Xj - i2wS) 


i2wJ3£ 


T e ‘ 1 ’ 2 


x,r ll>2 c. 


X l F 


CiC 2 


s' A 


( 1)2 


2w 2 J3 2 U ~ X t - i2wB) i2uSf? 
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r (l,2 c 2 

+ 2 a 2 r 2 " Cl 


2u 2 A 2 f 2 £ 


'Cl 

1 

T ‘ 1 » 2 L 

X 2 

f C 1 C 2 1 

1 4 - — 

CiC 2 

i T 

.2 

T s > 

J. s Jj 

. £ 

X T 

i2ufi, 

+ C ^ L 

i2wfif£ 


' 3T (2> \ 


c,c 


3 c 2 


1 -c 2 r 
+ - • - 


3y ) + i2wfl£ 
£.( 1)2 


' 1 ) 

- X 2 + ~ 


1 

+ - 


(\ 2 - X t )F 


£ (n t 2 - X t - i2uS) 


£ 2 w 2 fl 2 £ 

A , 

< 3T 


V 



X 2 

f C lC 2 ] 

CiC 2 

X 2 “ - 

~ C 2 L 


— 

1 + 

- + c : L 

i £ > 

„ 


. * 

i2uJ3y 

i2wBf£ 


- 1 


A .( 2 ) 


3y J 


c->c 


3 C 2 


i2 u& £ 


~ X 2 + — 


1 
r J 


+ - 
£ 


c, ' 


x 2 c 4 + c 5 


* c i L 



X 2 

r ^ i c 2 

Cl C 2 

1 ^ T 


. £ 

J- + 

i-2ufl j 

» C ^ Li 

i2«Af£ 


' 3T <2) X 


3y 


+ G, 


X ■> ( 0^2 \ 


[? 


1 + 


i2 uR j 


c.c 


l c 2 


12 u&f £ 


+ c ! L 


- 1 


(C. 30) 


where 


c, » c, 


1 ) 
T e 


.( 1)2 


c,r 


( 1)2 


c„ = - 


C c = 


2u 2 A 2 r 

X 1 ( X 2 “ \ i) 

( r X 1 2 - X t - i2ufl) 
(x 2 - X x ) F 
(rx 2 2 - X l - i2wfl) 


.( 1 ) 


( 1 ) 


c,r 


A ( 1 )x 


i uR j 


and 


G, = 


C 3 C 2 

f 1 

-• \ 4. __ 

1 

X — 

f r* \ 

c 4 

i2uB£ 

A 2 • 

x f ) 

T* — 
£ 

X 2 C 4 + c 5 

f > 


"* CoL 


are used. Combining Eqs. (C.15) and (C.30) yields Eg. (3.37) 
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Note that for the time dependent condition, the term T 8 tl>2 has 
to be substituted by half of its value as indicated in Eq. 

(3.13) . 

If the second order time independent solution is required, 
then the boundary equation is different from the previous one and 
it is given as follows: From Eq. (3.35), 


3T. 

£( 0 > + £ ( 2 ) 1 


9 t 

9y 


( o > 


- r 


ay 


a 2 r 


ay' 


= -r 


( i > 


3T, 


( l ) 


ay 


(C. 31) 


The characteristic equation in this case is given by 


x (1 - n) - 0 


(C. 32) 


Therefore, the homogeneous solution is of the form 

r h = D 0 + D x eY /f (C . 33 ) 

and the particular solution is 

r p = D 2 e x i Y + D 3 ye y/ ^ (C.34) 

Substituting Eq. (C.34) into Eq. (C.31) and investigating the 
coefficients gives the following relationships, i.e., 


D 2 

d 3 



.( 1 ) 2 , 


iwSf 


(C . 35 ) 


Now, the solution r is of the form 
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T ■ T h + T p * D 0 + Die y/r + 


£ x x ■* x i 


: X 1 y 


+ c< 


.< 2 ) 


A 

r' 


ye 


y/f 


i w sr j 

Since r = r 8 at y = 0 and r = 0 at y = -®, then 


D 0 - 0 


Dl = T 8 <2> - 


Finally, 


n x 2 * X x 


r * 


( 2 > 


+ c, 


.( 2 ) _ 


n x 2 - j 

A 

r 


e 


Y/f 


fX x 2 - X x 


: X 1 y 


A . ( 1 ) 2 \ 


iwBi” 


ye 


y/f 
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a r 

H 1 

( A f 

m < 2 > 

y 

eY /f h 


ay 

f 

i S £X x 2 - X x J 

£ x x 

+ 

C, 

f A ( 1 ) 2 X 

r‘ 2 > - 

e y/i> 

+ ° 2 

f 

A ( 2 ) 


X , F 


e i 


X lY 


iuftf 


.( 1)2 


iuflf 


ye 


y/f 


At y ■ 0 , 


9r 

3y 


y = 0 


T < 2 > - 


£ i i * — / 


X i F 


n x 2 - X x 


+ c. 


.( 2 ) 


and then from Eq. (3.10) 


(C. 36) 


(C. 37) 


(C.38) 


(C. 39) 


r * 1 ’ 2 ) 

iwSf i 
(C. 40) 
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Note that in actual calculation the term T s <1)2 has to be 
substituted by (1/2) T s ‘ 1 ’ ^Tg ' 1 ’ as indicated in Eq. (3.13). 
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2. Temperature Conditions at Flame Edge 

The temperature at the flame edge can be obtained by 
assuming that the flow very close to the flame is isentropic (see 
[21]). Writing this expression in terms of temperature and 
pressure yields 


p 

f p } 

7 

' T 

( 7/7- 1 ) 

p( 0 ) 

. ( 0 ) 
p y 


[ T (0) J 



(C. 43) 


Therefore, 


/ p \ ( 7- 1 ) /7 


,( 0 ) 


l P 


( 0 ) 


or 


|n 


,( o ) 


7-1 


jen 


l p 


( o ) 


(C. 44) 


Perturbing both sides yields 


|n 


1 + ( 


£ < i > 


<( o > 


+ e 


A 2 ) 


.< 0 ) 


7-1 


|n 


1 + e 


,< i > 


A 0 ) 


,( 2 > \ 


+ € ' 


,( 0 ) 


(C. 45) 

Since J2n(l + X) - x - (1/2) X 2 + (1/3) X 3 «... for |x| <1, Eq. 

(C. 45) can be rewritten as follows: 


,( i ) 


t 2 ) 


! (j( l ) 2 


• t 

T ( 0 ) 

l T (0) 

2 

T < 0 > 2 


+ t 2 

r p< 2 > 

1 

pill! 

- 

T C 

[ P <0> 

CM 

p<0)2 



+ 0 (€ 3 ) = 


+ 0 ( e 3 ) 


7-1 


.( 1 > 


,( 0 ) 


(C. 46) 
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Thus, for the first order term, we have 

iji(I) <y _ ^ p( 1 ) 

, (C.47) 

IJl(O) y p ( 0 ) 

and for the second order term 


$< 2 ) 

1 

£< 1 ) 2 7-1 

( p ( 2 > 

i 

1 H 

p< 1 ) 2 

T < 0 ) 

2 

IJ ( 0 ) 2 -y 

[ P' 0) 

2 

p(0)2 


Since DS/Dt =» 0, where S denotes the entropy, which means 
that the entropy is conserved for each fluid particle convected 
away from the flame, Eq. (C.47) is of the form 


[T tl, (y + Ay, x + ax, t) - T* 1 ) (y, x, 

( o > 

7-11 

= [p ( 1 1 (y + Ay, x + ax, t) 

7 P* 0 ’ 

Assuming that dT <1, / 9 y » aT (1) /9x and P (1> 
(flame edge), then 


t - At)] 

- P n, (y, x, t - At)] 

(C.49) 

is constant at y = | 


1 

T < o ) 


9T (1> 

+ i«T ll> 

ay 



i w 


p< i ) 


p( o ) 


Substituting Eq. (C.47) into Eq. (C.48) shows 


£( 2 > 7-1 

/ p< 2 ) 

1 

p< X ) 

2 

rp ( 0 ) -y 

l P< 0) 

27 

p( 0 ) 

2 

/ 


(C. 50) 


(C. 51) 


and using the same argument about the entropy conservation, 
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i( 0 ) 


3T 


( 2 > 


+ i2uT <2) 


3y 


7-1 


,< 2 ) 


i2 u 


pt 0 ) -yp< 0 ) 2 


iu ' 

pt 1 ) 2 


(C. 52) 


Combining Eqs. (C.50) and (C.52) yields 


1 (3T (I) 


,( o ) 


+ iluT 


3y 


— ilw 


7 - 1 P l 1 ’ 


, 10 ) 


(7 - 1) iw 


,( 1)2 


,(0)2 


or 


. 3T * 1 * 7 - 1 A (7 - 1) 

ilwT* 1 * + ilu P CI > iuP ( 

3y - 7 7 2 


(C. 53) 


which is equivalent to Eq. (3.39). 


3 • Fuel Species Condition at Interface 

Since the mass flux fraction of the fuel and oxidizer is 
fixed by the propellant composition, the species boundary 
condition at the gas-solid interface can be derived from Eq. 
(3.4). After neglecting the unsteady and reaction terms, the 
governing equation yields 


PU i Y /i - Y #ii - 0 (C. 54) 

For one-dimensional consideration, Eq. (C.54) can be integrated 
with respect to y and the results show 

1 dY 

m f = Y (C. 55) 

m dy 

where m is the mass flux and m f is the fuel mass flux fraction. 
Perturbing Eq. (C.55) yields 
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01 f 


( 0 ) 


+ fm f 
- [n 


<i> + € 2 m f ,2) = [Y 
A < 1 1 


( o > 


( 01 + €*"’ + 


+ (Y 
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Y ll) + e 2 Y l2) ] 
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-( o ) _ 
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£ 0 ) 1 


01 
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m (1> dY (0) 


.< o » 


01 '“' dy 
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+ e 


+ € 


y ( 2 > _ 


f m ( 1 1 2 

m ,2 > 

dY (0) y 

fll ,0 > 2 

_ t o > 

m J 

dy J . 


Y * l 1 _ 

1 j 

dY 1 1 > 

m<°> l 

dy 

1 

( dY (2> 

m (1) dY (1> 

m <0) 

i dy 

oi (0) dy 

1 

0(£ 3 ) 

(C. 56) 

01 ,0> 


Therefore, 
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( 0 > 


Ol <0) 

dy 
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m (1> 
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01 

dy 

0l (0) 

dy 

1 

r d ^<2> 
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dY <o> 

_( 0 1 
01 

l dy 

oi <0> 

dy 

1 m (1> dY * 1 * 

m ( 1 1 2 

dY t0> 

l 01 (0> 

dy 

01 <0)2 
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Since oi f 

Y < i > _ 

Y < 2 > _ 


(1) =» m ‘ 2 > 


(C. 57) 


0, and using Eqs. (B.12) - (B.13), we have 


' dY 1 1 ’ A dY 

. ^(01— m 111 

P s r C 1 1 S 


( 0 ) 


dy 


dy 


= 0 


(C . 58 ) 


rdY ,2> l 

' Cj 1 

\ dY (0 > 1 

_( 0 ) _ J m t 2 > . 


m (D2l 

P s r C 1 T s + 
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(C. 59) 
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Knowing that p s r l0> = p 0 v (0> = we have 


£< 1 > 


( 2 > _ 


dY (1> 


dy 

dY ( 2 * 


dY 


( o ) 


- c. 


- C. 


dy 

dY ,0 > 
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dy 


( i > 


T s l2> + 


(C. 60) 


< 1)2 


“ c x T s 


( l ) 


dY' 


+ T, 


( 1)2 


dY 


< o ) 


(C. 61) 


dy dy 

Combining the two equations yields Eq. (3.43). 

4. Normal Velocity Condition at Interface 

The boundary condition of the normal velocity component is 
decided using the mass balance condition at the interface. For 
the first order mass balance, we have 


.< l ) _ 


p‘0»v‘ I* + A p ( i>v ,0) 


Ps? 

Therefore, 

v* 11 - T (0> ( Ps r (1) - p l 1 , r ( 0) ) 
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(C. 62 ) 
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(C. 63) 


where 
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A 
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( l ) 
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and 

£( i) . « o > t i ) + ^ c i > T < o > 

are used. For the second order condition, 
Ps r « 2 > = p (0, v <2> + + £ ,2, v <0) 
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(C. 64) 


Combining Eqs. (C.63) and (C.64) yields Eg. (3.44). 


5. Parallel Velocity Conditions at Interface 

Since the parallel component of the velocity has to be zero 
at the actual solid surface in order to satisfy the no-slip 
condition, we can derive the higher order condition using the 
Taylor series expansion about the origin such that 
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3u + ‘ 0> 1 3 2 U + ‘ 0 ’ 

u + - u + <0> + 5y + ay 2 + ... (C.65) 

3y 2 3y 2 


Differentiating Eq. (C.65) with respect to t and multiplying p 
yields 


3u + 3u + <0) 1 3 2 u + 1 0 * 

p =» pv + + — /»v + $y+... (C. 66) 

3t 3y 2 3y 2 


For the first order condition, 
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( x > 
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( o > 


iwfl 3y 

For the second order condition, 
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(C. 67) 


Ill 


fl 8u (1) A B 3u (1> 

a £<2> + £< 1 > 

i2u 3y + l2w dy + 

fl 2 3 2 u (0) 

8 w 2 3y 2 

Combining Eqs. (C.67) and (C.68) yields Eq. (3.45). 


£( 1)2 

+ 


(C.68) 
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APPENDIX D 

DERIVATIONS OF FINITE ELEMENT EQUATIONS 


1. Zeroth Order Analysis 

The formulation of the zeroth order governing 
equations (3.16) is 

(ZA afl + ZB afi )T^ 0> = ZG a 

where 


ZA 


a B 


'ft 




ZB 


a fi 


^a, i ^ji, i 


'ft 


ZG a = 


T< ° ,n i dr + h 

/I 


ft 


( 0 ) 


(D.l) 


and 


w 


( o ) 
fl 


1 - T ( 0 * ) 2 


— Bz 
h 2 


,-E/T 


( 0 ) 


I ( 0 ) 


J fl 


2. Higher Order Analysis 

The finite element analogs for the higher order system 
are summarized below. 

Continuity 


[(iIu)FA a £ + FB a g + FC a g] Pq 


+ [FD 
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+ FE 
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( i ) 
B 
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( I > 


+ [FF afl + FH afl ]v fl = FG 0 


(D. 2 ) 


where 
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Momentum Equation . f X-component ^ 


[XA^fllPs 11 + C<iI“)XB a|5 + XC afl + XD afl + XE afl 


+ XF ia ]u'" + [XH„ + XI a .]vi>> + XJ aa p”’ - XG, 


tt B J B 


a B ' 


a B B 


'a 


(D.3) 



114 


where 


XA 




a 


( 0 ) ( 0 ) 
u 3i dA 


XB 


a B 


A 


< 0 ) 

taLt&tyPy dA 


XC 


a B 


( o ) ( o ) 


*a*A, j *y*sPy u 3j dA 


A 


XD, 


a fl 


( 0 ) ( 0 ) 

*,j Py u ji dA / (j = 1) 


A 


XE afl - Pr 


SI 


dfl 


XF 


a B 


1 

- Pr 
3 


SI 


XH 


a fl 


a 


^a, i ^A/ j dfl / ( j 88 1) 


( 0 ) ( 0 ) 

, j P >7 u 3i dfl , (j -■ 2) 


A 


XI - - Pr 
3 


SI 


dA , (j — 2) 


XJ 


0,6 7M b 2 d 


A 


i dA 


XG„ 


A 


( 0 ) A ( 1 ) ( 1 ) 
U 7j u 3i dA 


A 


A ( 1 ) ( 0 ) ( 1 ) 

, i Pr u 7 j u s i dA 


A 


A ( 1 ) A ( 1 ) (0) 

*a* B Ma, j PR u 7j u si dA + Pri 


A ( I > * 

u i,j n j*a dr 



115 


1 
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u j,j n i*a <*r 
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Here, all i = 1. 

Momentum Equation. ( Y-comoonent^ 
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Energy Equation 
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Species Equation 

( i ) ( i ) ( I ) ( i > 

[SA. fl + SB a# ], s + S C<lfl U g + SD ag v g + SE ag T g 


( I ) 


+ [ (il«) SF aj j + SH aft + SI aft + SJ aft ]Y fl = SG 




a B 


a B J "B 


where 


SA 


a B 


( o ) ( o ) 


*a*fl*Y*» , i u Yi Y J dn 


SB 


a B 


d A 


) 

f 2w <0> ' 


( 0 ) 

A 

p 


( o ) 

T s dA 


(D. 6) 



118 


SC, 


( 0 ) ( 0 ) 




A 


tat&tyti , i Py Y « dft , (i - 1) 


SD 


a fl 


A 


( 0 ) ( 0 ) 

fatbits, i Py Y 8 / (i - 2) 


SE 


a fl 


'XI 


*a* fi *7 


r Ew (0) 


l T 


( 0)2 


dA 


SF 


a fl 


'A 


( o ) 

^ a.^ fo^y Py di ^ 


SH 


( 0 ) ( 0 > 


a fl 


/ i*y*lPy Ug { dA 


A 


SI 






ft 


SJ 


a fl 


A 


*a*S*7 


t 2w (0) 


l Y 


( o ) 


dA 


SG a = - 


A 


( 0 ) ( 1 ) ( l ) 

,iPfl U 7i Y S d ^ 


A 


( 1 ) ( 0 ) ( 1 ) 

/ i Pjj u 7i ^8 d ^ 


a 


A 


( 1 ) ( 1 ) ( 0 ) 
i P q u 7 i ^8 d ^ 


'A 


A 




7*8 


Ew 


( 0 ) 


,(0)3 


( E \ 

- 1 


2T 


( 0 ) 


( 1 ) ( 1 ) 

T j dA 


fl 


f 2Ew 


( o ) 


< 0 > T ( 0 ) 2 


( 1 ) ( 1 ) 

Py Tg dA 


6 



119 



n 


i 

" a 

— lu 


( 2Ew 


( 0 ) 


** *Q*y*s 


[ y ‘ 0 > t 1 0 * 2 j 

' 4w ( 0) 


Cl ) ( 1 ) 


dfi 


fl 


V p' 0) Y C0) 


/ w c o ) ^ 


Cl) Cl) 

Py Yj dfl 


fl 


v p 


( 0)2 


< 1 ) Cl ) 

Py p S dfl 


f W (0) A 


( 0)2 


Cl ) ( 1 ) 


dfl 


fl 




*°-*r *ypR 


c i ) 


c i ) 

Yy dfl 


